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We examine five diflterent popular rigid water models (SPC, SPCE, TIP3P, TIP4P and TIP5P) 
using molecular dynamics simulations in order to investigate the hydrophobic hydration and inter- 
action of apolar Lennard- Jones solutes as a function of temperature in the range between 275 K 
and 375 K along the 0.1 MPa isobar. For all investigated models and state points we calculate the 
excess chemical potential for the noble gases and Methane employing the Widom particle insertion 
technique. All water models exhibit too small hydration entropies, but show a clear hierarchy. 
TIP3P shows poorest agreement with experiment whereas TIP5P is closest to the experimental 
data at lower temperatures and SPCE is closest at higher temperatures. As a first approximation, 
this behaviour can be rationalised as a temperature-shift with respect to the solvation behaviour 
found in real water. A rescaling procedure inspired by information theory model of Hummer et al. 
(Chem.Phys. 258, 349-370 (2000)) suggests that the different solubility curves for the different mod- 
els and real water can be largely explained on the basis of the different density curves at constant 
pressure. In addition, the models that give a good representation of the water structure at ambient 
conditions (TIP5P, SPCE and TIP4P) show considerably better agreement with the experimental 
data than the ones which exhibit less structured 0-0 correlation functions (SPC and TIP3P). In the 
second part of the paper we calculate the hydrophobic interaction between Xenon particles directly 
from a series of 60 ns simulation runs. We find that the temperature dependence of the association 
is to a large extent related to the strength of the solvation entropy. Nevertheless, differences between 
the models seem to require a more detailed molecular picture. The TIP5P model shows by far the 
strongest temperature dependence. The suggested density-rescaling is also applied to the chemical 
potential in the Xenon-Xenon contact-pair configuration, indicating the presence of a temperature 
where the hydrophobic interaction turns into purely repulsive. The predicted association for Xenon 
in real water suggest the presence a strong variation with temperature, comparable to the behaviour 
found for TIP5P water. Comparing different water models and experimental data we conclude that 
a proper description of density effects is an important requirement for a water model to account 
correctly for the correct description of the hydrophobic effects. A water model exhibiting a density 
maximum at the correct temperature is desirable. 

PACS numbers: 02.70.Ns,61.20. Ja,61.25.Em,82.60.Lf 



I. INTRODUCTION 



Nonpolar solutes show a strong tendency to aggre- 
gate when dissolved in water due to the relatively strong 
water- water interaction in comparison to the weak solute- 
water interaction However, in addition to such 
energetical considerations the hydration entropy of small 
simple solutes is found to be negative, which is usually 
explained as increased ordering of the molecules in the 
hydration shell 0, . This is the characteristic feature 
of the so called hydrophobic hydration of small apolar par- 
ticles TJ. As a consequence, with increasing tempera- 
ture the association of hydrophobic particles is found to 
be enhanced in order to minimise the solvation entropy 
penalty H, 0] . This entropy-driven association process 
is usually referred to as hydrophobic interaction and has 
been subject of numerous simulation studi es l3.I^ITollTlL 
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Hydrophobic effects are considered to play an impor- 
tant role concerning prot ein stability and other self as- 
sembly phenomena 4, 28, "295. The strength of hydropho- 
bic interactions, however, might as well determine largely 
; the temperature-range where the protein remains physio- 
logically functional before thermal unfolding takes place. 
^ In addition, the weakening of hydrophobic interactions 
[ at lower temperatures is presumably of importance for 
the understanding of the cold denaturation of proteins 
! [30 . 31, 32]. This is probably of particular importance 
; for systems that show an entropy driven configurational 
"folding" transition such as the "hydrophobic" collapse 
of polymeric elastin [33L l3^ and small tropoclastin-like 
I oligo-peptides around 300 K — 330 K. The present 
status of conceptual understanding of hydrophobic hy- 
j dration and hydrophobic interaction has been reviewed 
; recently by Pratt and Southall et al. [s^], whereas 
Smith and Haymet 38] give a tutorial overview over the 
currently available computational methods. 

Recent methodological improvements in simulation 
techniques, such as the "parallel tempering" approach 
'7de/''tpaii5/, replica exchange molecular dynamics") give rise to 
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the possibility of explicitly studying the thermal fold- 
ing/unfolding transition of solvated proteins [s^ l40l| . 
Molecular dynamics simulations of proteins in aqueous 
solution, however, still largely depend on the use of sim- 
ple model potentials for water. Therefore it might be in- 
teresting to compare different popular rigid water models 
with respect to hydrophobic interactions and to show how 
their strengths might be related to the behaviour of other 
known thermodynamical and structural properties of the 
model liquids. For this purpose we employ the three cen- 
ter point charge models proposed by the Berendsen group 
(SPC 113 and SPCE 42]) as weU as three candidates 
representing the TIP-family of water models according 
to Jorgensens group: The three center TIP3P and four 
center TIP4P 43\, as well as the recently proposed five 
center TIP5P model 

The most simple and as well most often studied hy- 
drophobic model system in this context is, of course, hy- 
drophobic Lennard- Jones particles dissolved in water |45l |. 
The outline of the paper is therefore the following. 

First we would like to examine the the performance of 
the five different models with respect to the hydrophobic 
hydration behaviour as a function of temperature with 
a focus on the physiologically important temperature 
range between 275 K and 375 K. This is done a in the 
spirit of the paper by Guillot and Guissani |46j | where we 
calculate the chemical potentials and solvation entropies 
for Lennard Jones particles representing the noble gases 
and Methane and compare them with experimental data. 
Simulation runs of 20 ns allow an accurate determination 
of the excess chemical potential using the Widom particle 
insertion method. 

In the second part of the paper we study the associa- 
tion behaviour of hydrophobic particles (only Xenon) as 
a function of temperature for all five water models by 
conducting long (60 ns) MD simulation runs of solutions 
of the hydrophobic particles in water, as suggested re- 
cently by the work of Ghosh et al. |2a . l27| . The potential 
of mean force (PMF) between the hydrophobic particles 
can then be obtained directly from the pair distribution 
functions. Thus we are able to calculate the hydropho- 
bic interaction for the different models as a function of 
temperature. 

A very recent study on lattice models by Widom et al. 
[47| suggests a linear relationship between the strength 
of the hydrophobic hydration and and the hydrophobic 
pair interaction. Therefore it might also be interesting to 
provide accurate data for both on a set of realistic water 
models. 



II. METHODS 

A. MD Simulation details 

We employ molecular dynamics (MD) simulations in 
the NPT ensemble using the Nose-Hoover thermostat 
[49I [soi l and the Rahman-Parinello barostat [sil [H^l] with 



coupling times tt = 1.5 ps and Tp = 2.5 ps (assuming the 
isothermal compressibility to be xt = 4.5 lO^^bar"^), 
respectively. The electrostatic interactions are treated 
in the "full potential" approach by the smooth particle 
mesh Ewald summation | 53i] with a real space cutoff of 
0.9 nm and a mesh spacing of approximately 0.12 nm and 
4th order interpolation. The Ewald convergence factor a 
was set to 3.38 nm^^ (corresponding to a relative accu- 
racy of the Ewald sum of 10^^). A 2.0 fs timestep was 
used for all simulations and the constraints were solved 
using the SETTLE procedure AH simulations re- 

ported here were carried out using the GROMACS 3.1 
program |55l Is^ . Statistical errors in the analysis were 
computed using the method of Flyvbjerg and Petersen 
[s?! . For all reported systems and different statepoints 
initial equilibration runs of 1 ns length were performed 
using the Berendsen weak coupling scheme for pressure 
and temperature control TT = Tp = 0.5ps [s^. 



In order to determine the excess chemical potential of 
the hydrophobic particles we performed a series of sim- 
ulations generally using 256 water molecules for all five 
different water models SPC, SPCE, TIP3P, TIP4P and 
TIP5P (model parameters are given in Tabled). How- 
ever, the excess chemical potential is known to be sensi- 
tive to finite size effects [53 ■ In order to estimate this in- 
fluence, we additionally carried out simulations contain- 
ing 500 and 864 water molecules but only for the SPCE 
model. All model systems were simulated at five differ- 
ent temperatures 275 K, 300 K, 325 K, 350 K and 375 K 
at a pressure of 0.1 MPa. Each of these simulations ex- 
tended to 20 ns and 2x10"* configurations were stored for 
further analysis. To determine the hydrophobic interac- 
tion between Xenon particles (for the model parameters 
see Tabled we use MD simulations containing 500 water 
molecules and 8 Xenon particles employing the same sim- 
ulation parameters outlined above. Again, the five dif- 
ferent water model systems are studied at 275 K, 300 K, 
325 K, 350 K and 375 K at a pressure of 0.1 MPa. Here, 
runs over 60 ns were conducted, while storing 7.5 x 10^ 
configurations for further analysis. The simulations pro- 
tocols showing the obtained densities and average poten- 
tial energies are given in Table HTl 



Concerning the model parameters we would like to em- 
phasise that the parameters describing the noble gases 
(see Table|2for details) used here were fitted to reproduce 
their pure component properties. Hence a perfect match- 
ing of the solubilities with the experimental data cannot 
be expected. However, since we are interested in compar- 
ing different water models, taking these parameters is the 
preferred procedure since they should work for all models 
equally good (or bad). The water/gas cross terms were 
obtained applying the standard Lorentz-Berthelot mixing 
rules according to (jy = {an + Ujj) /2 and = ^^a^jj- 
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Model 


tr/A 




g/e 


SPC 


3.1656(0) 


78.2(0) 


0.41(H) 


SPCE 


3.1656(0) 


78.2(0) 


0.4238(H) 


TIP3P 


3.1506(0) 


76.58(0) 


0.417(H) 


TIP4P 


3.1536(0) 


78.08(0) 


0.52(H) 


TIP5P 


3.12(0) 


80.56(0) 


0.241(H) 


Ne 


3.035 


18.6 





Ar 


3.415 


125.0 





Kr 


3.675 


169.0 





Xe 


3.975 


214.7 





CH4 


3.730 


147.5 






TABLE I: Lennard-Jones potential parameters and partial charges describing the water-water and solute-solute pair interactions. The 



solute-water cross parameters are deduced from the Lorentz-Berthelot mixing rule s: a - . 
information on t he geom etry of the water models we refer to the original references |41L |4 
taken from Refs. El El. 



,-)/2, 



For further 



The solute-solute-parameters were 



B. Infinite dilution properties 

Usually the solubility of a solute is measured by the 
Ostwald coefficient i'^^ =Pb/Pb' where and are 
the number densities of the solute in the liquid and the 
gas phase, respectively, when both phases are in equilib- 
rium. Here A denotes the solvent and B indicates the 
solute. Equilibrium between both phases leads to a new 
expression for L'/^, namely, 



exp 



(1) 



where (3 = 1/ kT and g and /i^^ ^ denote the excess 
chemical potentials of the solute in the liquid and the gas 
phase, respectively. When the gas phase has a sufficiently 
low density, then s ~ 0, hence L'/^ becomes identi- 
cal to the solubility parameter 7^ = exp [— /3^ea; b] ■ ^or 
our study, covering the temperature range between 275 K 
and 375 K, the excess chemical potential of apolar solutes 
in the gas phase can be practically considered to be zero 
(see Table 3 in Ref 0| ) . The chemical potential of a so- 
lute can be obtained from a constant pressure simulation 
(NPT-Ensemble) of the pure solvent using the Widom 
particle insertion method ,60, 61] according to 



Mb 



In 



A3 



-/3-1 In iLL^^^^+i ^''P^-^ 



(V) 



(2) 



where AU = t/(s^+^ L) - U{s^;L) is the potential en- 
ergy of a randomly inserted solute {N + 1)- particle into 
a configuration containing N solvent molecules. The 



Si = L^^ri (with L = V^^^ being the length of a hy- 
pothetical cubic box) are the scaled coordinates of the 
particle positions and / sn+i denotes an integration over 
the whole space. The brackets (. . .) indicate isothermal- 
isobaric averaging over the configuration space of the A'^- 
particle system (the solvent). A represents the thermal 
wavelength of the solute particle. The first term /j,'^ g is 
the ideal gas contribution of the solute chemical poten- 
tial at an average solute number density (p^) = 1/ (V^) 
at the statepoint (T, P) . We would like to point out that 
the definition of the in Eq. |2] assumes that solute and 
solvent particles are of different type and hence distin- 
guishable. Since we are considering water at relatively 
low temperatures, the volume fluctuations are compara- 
bly small. Hence the obtained values for /iea; = /Xg^. ^ are 
practically identical to the values obtained from constant 
volume simulations at the same statepoints with differ- 
ences due to the fluctuating volume < 0.02kJmol~^ fore 
the temperature range considered in our study. The en- 
tropic and enthalpic contributions to the excess chemical 
potential can be obtained straightforwardly as tempera- 
ture derivative according to 



dT 



and 



Me 



T Se 



(3) 



and the isobaric heat capacity contribution according to 



-T 



QT2 



(4) 



As an alternative to the Ostwald coefficient, the sol- 
ubility of gases is often expressed in terms of Henry's 
constant kn- The relationship between Henry's constant 
and the solubility parameter 7^ in the liquid phase is 
given by ||62] 



kn = p!4^T/7^ 



(5) 
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Model 


T/K 


(p) /kgm--' 


(E) /kJ mor^ 


(p) /kgm-3 


(E) /k,]mori 


SPC 


275 


993.3 ±0.2 


-43.005 ± 0.003 


1070.0 ±0.1 


-42.470 ± 0.001 




300 


977.1 ± 0.1 


A -1 r o /I 1 n nnn 

—41.534 ± 0.002 


1 n r n n i n 1 

1050.2 ± 0.1 


A r\ ni — 1 1 n nni 

—40.951 ± 0.001 




325 


956.7 ± 0.1 


A r\ nn/^ i n nni 

—40.096 ± 0.001 


1 nn/^ f 1 n 1 

1026.5 ± 0.1 


on A 1 n nni 

—39.472 ± 0.001 




350 


933.8 ± 0.1 


oo /^n^ 1 n nnn 

—38.687 ± 0.002 


nnn r' i n t 

999.5 ± 0.1 


oonin 1 nnni 

—38.019 ± 0.001 




375 


907.8 ± 0.1 


o^ 1 n nno 

—37.277 ± 0.003 


n n n i n t 

969.2 ± 0.1 


o r^^i 1 n nni 

—36.571 ± 0.001 


SPCE 


275 


1009.5 ±0.2 


-48.148 ±0.004 


1089.9 ±0.1 


-47.594 ± 0.003 




300 


998.6 ± 0.2 


J r' ^n 1 n nnn 

—46.572 ± 0.002 


1076.1 ± 0.1 


/ir'n/^o t nnni 

—45.968 ± 0.001 




325 


no o o 1 n 1 

983.8 ± 0.1 


-ir" dn T' 1 nnm 

—45.066 ± 0.001 


1057.5 ± 0.1 


A A /<no 1 n nnn 

—44.408 ± 0.002 




350 


965.3 ± 0.2 


.1 o r n o 1 n nnn 

—43.598 ± 0.002 


1 no r 1^ t n 1 

1035.7 ± 0.1 


/in on'7 1 n nni 

—42.897 ± 0.001 




375 


r\ A A 1 n n 

944.3 ± 0.2 


A c\ -x v A 1 n nnn 

—42.154 ± 0.002 


1 n 1 n n i n i 

1010.2 ± 0.1 


A 1 /I n f 1 n nnn 

—41.405 ± 0.002 


SPCE (500 Mol.) 


275 


1009.0 ±0.1 


-48.147 ±0.003 








300 


nno A \ r\ ^ 

998.4 ±0.1 


Ad c Tr' 1 c\ r\r\ 1 

—46.576 ± 0.001 








325 


no o o 1 n 1 

983.3 ± 0.1 


A r" n/^n i n nni 

—45.062 ± 0.001 








350 


n J n 1 n "1 

964.9 ± 0.1 


/I o r'on 1 n nn-i 

—43.589 ± 0.001 








375 


n .1 o o 1 n T 

943.8 ± 0.1 


A c\ '\ A <\ 1 n nnn 

—42.142 ± 0.002 






SPCE (864 Mol.) 


275 


1009.0 ±0.1 


-48.152 ±0.002 








300 


nno n 1 n 1 

998.2 ± 0.1 


/I f >~r'~r 1 n nni 

—Ab.bi" ± 0.001 








325 


n o o o 1 n "1 

983.3 ± 0.1 


/ir^ n/^/"* 1 n nnn 

—45.066 ± 0.002 








350 


n/^ A a 1 n 1 

964.8 ± 0.1 


/I o fnn 1 n nni 

—43.592 ± 0.001 








375 


943.6 ± 0.1 


/in 1 /in 1 n nni 

—42.149 ± 0.001 






TIP3P 


275 


1005.2 ±0.1 


-41.300 ±0.002 


1079.8 ±0.1 


-40.753 ± 0.001 




300 


no A r' 1 n 1 

984.6 ±0.1 


on nni i n nni 

—39.921 ± 0.001 


1 n r' r' o t n n 

1055.8 ± 0.0 


on onn ) n nni 

—39.329 ± 0.001 




325 


n/^n n 1 n 1 

960.9 ± 0.1 


OO f rr" 1 n nnn 

—38.565 ± 0.002 


1 nn o /I inn 

1028.4 ± 0.0 


o no*? 1 n nni 

—37.937 ± 0.001 




350 


no ,1 z"' 1 n n 

934.6 ± 0.2 


o nnn i n nnn 

—37.229 ± 0.002 


nn^ 1^ 1 n 1 

997.7 ± 0.1 


o r'/^n 1 n nni 

—36.560 ± 0.001 




375 


nnr" a inn 

905.4 ± 0.2 


o f o o o 1 n nnn 

—35.888 ± 0.002 


n o n 1 n 1 

963.9 ± 0.1 


or' ion 1 nnni 

—35.182 ± 0.001 


TIP4P 


275 


1005.3 ±0.1 


-42.839 ±0.003 


1089.0 ±0.1 


-42.957 ±0.002 




300 


nno r" 1 n 1 

993.5 ± 0.1 


A -1 n o 'T' 1 n nno 

—41.237 ± 0.003 


1 n'^r' n 1 n 1 

1075.2 ± 0.1 


/II c\ A r\ 1 nnni 

—41.249 ± 0.001 




325 


976.5 ± 0.1 


on/^'ni^ 1 nnnn 

—39.697 ± 0.002 


1 n r' r' n i n i 

1055.0 ± 0.1 


on/^nj 1 nnni 

—39.624 ± 0.001 




350 


955.1 ± 0.1 


OO nn/^ 1 n nnn 

—38.206 ± 0.002 


1 nnn n f n n 

1029.9 ± 0.0 


onnr^i t nnni 

—38.051 ± 0.001 




375 


nnn a i n i 

929.4 ±0.1 


o 'Tc^r' 1 n nno 

—36.726 ± 0.003 


1 nnn f i n i 

1000.5 ± 0.1 


o r'no 1 n nni 

—36.508 ± 0.001 


TIP5P 


275 


987 8 ± fl 9 


—49 744 ± n 009 


1 OfiQ 5 ± 1 


—49 379 ± n04 




300 


982.6 ±0.2 


-40.132 ±0.004 


1058.6 ±0.1 


-39.650 ± 0.002 




325 


964.1 ±0.1 


-37.910 ±0.003 


1034.2 ± 0.1 


-37.314 ±0.002 




350 


936.6 ±0.1 


-35.887 ± 0.005 


998.5 ±0.1 


-35.248 ± 0.002 




375 


902.9 ±0.1 


-33.995 ± 0.002 


960.3 ±0.1 


-33.258 ± 0.002 



TABLE II: Average densities p and configurational energies E (per molecule) describing the examined statepoints at p = 0.1MPa. Left 
columns: pure water simulations. The system size is 256 molecules except for the SPCE simulations indicated. Each statepoint has been 
simulated for 20 ns. Right columns: The system size is 500 water molecules plus 8 Xenon atoms. All systems were simulated for 60 ns. 



where is the number density of the solvent. We use 
this relation in order to compare the experimental with 
the simulation data. 

The thermodynamic solvation properties discussed in 
this paper belong to the so called number density scale. 
In the experimental literature [6^ l6^ IbsI l6^. however, 



the properties are often discussed on the mole fraction 
scale with the solvation free energy being 

AG = -RT\iikH, (6) 

where Henry's constant is expressed in bars [g^, 113 ■ 
When comparing with experimental data, care must been 
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taken to which scale the discussed properties (solvation 
enthalpies, entropies and heat capacities) belong. Ther- 
modynamic properties determined on the mole fraction 
scale contain additional terms depending on the thermal 
expansivity of the liquid |6^ . 

In order to perform the calculation most efficiently 
we have made use of the excluded volume map (EVM) 
technique ,68. .6ft | by mapping the occupied volume 
onto a grid of approximately 0.2 A mesh- width. Dis- 
tances smaller than 0.7 x cTy with respect to any solute 
molecule (oxygen site) were neglected and and the term 
exp{—(3 AU) taken to be zero. With this setup the sys- 
tematic error was estimated to be less than 0.02 kJ mol~^. 
Although the construction of the excluded volume list 
needs a little additional computational effort, this sim- 
ple scheme improves the efficiency of the sampling by al- 
most two orders of magnitude. For the calculation of the 
Lcnnard-Jones insertion energies AU we have used cut- 
off distances of 9 A in combination with a proper cut-off 
correction. Each configuration has been probed by 10'^ 
successful insertions (i.e. insertions into the free volume 
contributing non- vanishing Boltzmann- factors) . 

For the case of Xenon we would also like discuss the 
effect of having a polarisable solute. Therefore we cal- 
culate as well an additional polarisation term according 
to 

AU = AULj + AUpoi (7) 

with 

AC/po/ = ~la\F\^ , (8) 

where a = 4. 11 A is the Xenon polarisability and F is the 
electric field created by all water molecules at the location 
where the particle is inserted. F is evaluated using the 
classical Ewald summation technique with a Ewald con- 
vergence factor of 2.98 nm^^ (corresponding to a relative 
accuracy of the Ewald sum of « 10~*) in combination 
with a real space cut-off of 9 A and a reciprocal lattice 

cut-off of \kjnax\'^=25. 

We have tested our calculations by recalculating the 
chemical potential for various noble gases and Methane 
at exactly the same statepoints as reported by Guillot 
and Guissani while explicitly taking the solute polar- 
isability into account. For higher temperatures (> 473 K) 
we can quantitatively reproduce their data, whereas for 
the lower temperatures (the temperature range of our 
study) and larger particles (Methane, Xenon) differences 
occur, but which are qualitatively in accordance taking 
the the estimated error of their relatively short calcula- 
tions into account. 

When studying particularly dense liquids, the Widom 
methods is known to fail |6lj |. In order to confirm the 
applicability of the Widom method we have performed 
additional checks on the accuracy of the obtained chem- 
ical potentials by comparing with results according to 
the overlapping distribution method (see section ^ for 
details). 



C. Hydrophobic Interaction 

We use simulations containing 500 Water molecules 
and 8 Xenon particles to study the temperature depen- 
dence of the association behaviour of Xenon. The hy- 
drophobic interaction between the dissolved Xenon par- 
ticles is quantified in terms the hydrophobic cavity po- 
tential w{r) 0. The w{r) is obtained by inverting the 
Xenon- Xenon radial distribution functions g{r), to get 
the potential of mean force (PMF), and subtracting the 
Xenon-Xenon pair interaction potential 

wir) = ~kTlng{r) - Vxc-Xeir) . (9) 

We use temperature derivatives of quadratic fits of 
w{r, T) to calculate the enthalpic and entropic contri- 
butions to at each Xenon-Xenon separation r. For the 
fits aU five temperatures 275 K, 300 K, 325 K, 350 K and 
375 K were taken into account. The entropy and enthalpy 
contributions are then obtained as 

and 

h{r)^w{r)+Ts{r) . (11) 

In addition, the corresponding heat capacity change rel- 
ative to the bulk liquid is vailable according to 

III. RESULTS AND DISCUSSION 
A. Density-Curves 

The simulated statepoints are summarised in Table ITU 
Here the obtained average densities and potential en- 
ergies are given. In addition, cubic fits of the density 
with respect to temperature were performed and the ob- 
tained fitting-parameters, which will be used later for 
the evaluation of the hydrophobic solvation properties, 
are listed in Table IIIII The density-fits, as well as the 
original data are shown in Figure ^ A rather striking 
difference between the water models is course the loca- 
tion of the density maximum. The TIP5P model was 
actually parameterised to yield a density maximum at 
exactly the experimental temperature and density |44| . 
However, Mahoney and Jorgensen used a cutoff of 1.2 
nm for the water-water interaction without any correc- 
tions. Since we apply the Ewald technique summation 
here, we obtain densities which are consistently about 
2% smaller than the values reported by Mahoney and 
Jorgensen, which is, however, in accord with the obser- 
vations of Lisal et al. [t^ ■ 

We would like to point out that our density-fits, al- 
though obtained at much higher temperatures, reproduce 
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Model 



Po/kgm 



pi/kgm ^K" 



P2/kgm ^K" 



pa/kgm ^^K" 



SPCE 

SPC 

TIP3P 

TIP4P 

TIP5P 



0.65860 10^ 
0.78387 10^ 
0.99460 10^ 
0.66286 10^ 
-0.83369 10^ 



0.34537 10^ 
0.25154 10^ 
0.95254 10° 
0.32444 10^ 
0.16033 10^ 



-0.99533 10"^ 
-0.78904 10"^ 
-0.37325 10-2 
-0.86816 10-2 
-0.44692 10-^ 



0.73989 10-*= 
0.55036 10-^ 
0.14878 lO-'^ 
0.51335 10-5 
0.38097 lO-"* 



TABLE III: Polynomial fits of the densities obtained for the pure water simulation series at 0.1 MPa pressure: p(T) = po 

P3 r3. 



-P1T + P2 T2 



almost quantitatively the location of the density maxima 
of SPCE water [llllllli and TIP4P water The 
density maxima for SPC and TIP3P are perhaps shifted 
to even lower temperatures (44i . i76j and thus might even 
lie below a possible glass transition [t^ . 

We would like to emphasise that in line with the tem- 
perature dependence of the density, the different mod- 
els show the same hierarchy with respect to their ability 
to reproduce waters structural features as shown by the 
work of Head-Gordon et al. ^78, 79, 80] . At ambient 
conditions the TIP5P model agrees quantitatively with 
experimental the 0-0-structure factor, SPCE and TIP4P 
agree well, whereas SPC and TIP3P are completely lack- 
ing the second peak in the 0-0 pair correlation function 
[73. Isof . Hence both structure and density suggest con- 
sistently that TIP4P and SPCE might be considered to 
reflect states of water perhaps at slightly elevated temper- 
ature, whereas TIP3P and SPC correspond (structurally) 
to states of water at much higher temperature. It should 
be mentioned that all models discussed here show a sub- 
stantial disagreement with ex per iment with respect to 
the OH and HH correlations [Tg- This might be par- 
tially attributed to the fact that the models used here 
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FIG. 1: Density of water as a function of temperature for the 
different water models employed in the present study. The thick 
solid line represents experimental data according to Wagner and 
Prufi I71I . The thin lines represent cubic fits to the simulated data 
(see table ITTTl for the fitting-parameters). 



are rigid. But, since as well flexible and polarisable mod- 
els as well as recent Car-Parinello simulations show the 
same tendency [t^, the reason for this is at present not 
clear. 

Both, solvent structure and density are considered to 
be of importance for the hydrophobic effects. Since a 
clear hierarchy of the quality of the water models com- 
pared with real water is observed with respect to those 
properties, the hydrophobic hydration and interaction 
properties might as well be affected in such a systematic 
way. 



B. Hydrophobic Hydration 

The calculated excess chemical potentials for the noble 
gases and Methane are shown in Figure|21and are given in 
Table llVl The experimental data for Xenon and Methane 
shown in Figure |21 have been calculated using Henry's 
constants according to Fernandez-Prini and Crovetto |81| 
employing the water densities for the 0.1 MPa isobar re- 
ported by Wagner and Prufi [tJ . The thin lines in Figure 
121 represent fits of the data to the modified information 
theory (MIT) model, discussed in section UlI CI The cor- 
responding fit parameters are given in Table FVl 

The errors for the excess chemical potentials were cal- 
culated using the method of Flyvbjerg and Peterson [s^ 
and are indicated in Table IIVI A further check on the 
consistency and accuracy of the data has been performed 
by application of the method of overlapping distribution 
functions discussed in section^ We would like to point 
out that the pLf,x data for polarisable Xenon obtained 
here do only qualitatively agree with the data of Guil- 
lot and Guissani However, test-calculations (results 
not shown here) performed for selected statepoints given 
by Guillot and Guissani at higher temperatures jig do 
quantitatively agree, suggesting that the differences at 
lower temperatures are due to the relatively large statis- 
tical error in their calculations and have to be attributed 
to their rather short simulation runs. Moreover, the ob- 
tained value for the excess chemical potential of Methane 
in TIP4P-water at BOOK of (9.78 ± 0.1) kJmol"^ agrees 
well with the value of (9.79 ± 0.21) kJ mol"^ obtained by 
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Shimizu and Chan |23 for 298 K and 1 atm. 

The SPCE simulations containing larger numbers of 
molecules (500 and 864) indicate that the excess chem- 
ical potentials obtained from 256 water molecules runs 
are subject to a systematic error of ^ex due to finite size 
effects. The observed system-size dependence is presum- 
ably due to the restriction of volume fluctuations up to a 
certain maximum wave-length according to the presence 
of periodic boundary conditions. As a consequence the 
values from the 256 molecule simulations are lying sys- 
tematically too high. As shown in Tabic Hvl this contri- 
bution is relatively small for Neon (« O.lkJmoP^), but 
increases with particle size and leads to approximately 
O.SkJmoP^ and 0.5kJmon^ too high excess chemical 
potentials for Methane and Xenon, respectively. How- 
ever, this difference is found to be only weakly tempera- 
ture dependent. Therefore the proposition of our paper, 
a focus on the temperature dependence, is not affected. 
The shift with respect to the true excess chemical po- 
tential can be expected to be in the same range for all 
models 's^ since the compressibilities of the different wa- 
ter models have the same order of magnitude (between 
4.6x 10-^bar"^(SPCE) and 6.0 x f 0"^ bar^^(TIP3P) at 
300 K). 

Figure El reveals that in all cases the simulated excess 
chemical potentials behave qualitatively like the exper- 
imental data in the sense that the excess chemical po- 
tential is positive and increases with temperature indi- 
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FIG. 3: Entropy contribution to tlie excess chemical potential 
Sex for Methane and Xenon. The thick dashed lines represent th e 
experimental data according to Fernandez Prini and Crovetto l8ll . 
The cl ose d circles are according to the experimental data of Rettich 
et al. |64| . Th e employed water densities are according to Wagner 
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FIG. 2: Excess chemical potential of Methane and Xenon in water 
as a function of temperature. The ex perimental data (heavy dashed 
line) are according to Refs. I7ll l8ll . The thin lines represent fits 
to the simulated data employing the MIT-model. 



eating a negative entropy of solvation. Please note also 
that around 275 K the simulated excess chemical poten- 
tials for all models except for the TIP5P model are lying 
rather close together, being situated about 2.5kJmol~^ 
(Methane) and 3.5kJmol~^ (Xenon) above the experi- 
mental data. With increasing temperatures, the differ- 
ences between the values corresponding to the different 
water models start to increase, already suggesting dif- 
ferently strong solvation entropies. The overall change 
of /iea; with temperature is smallest for the TIP3P and 
SPC model, larger for TIP4P and SPCE and most ex- 
treme for the TIP5P model. Please note also that for the 
TIP3P and TIP5P models the curves for Methane and 
Xenon suggest the presence of a maximum of fi^^ close 
to 375 K, whereas the experimental maximum is located 
at much higher temperatures. 

In Figure 13 the entropy contribution to the excess 
chemical potentials are shown for the different models 
as well as for the experimental data. The thin lines 
were obtained as numerical derivatives of MIT-model fits 
whereas the symbols represent finite differences of the 
/.tex-data given in Table ITVl We would like to point out 
the the experimental data shown here belong to the num- 
ber density scale and are therefore smaller than the values 
given in the paper of Rettich et al. |64j . corresponding 
to the mole fraction scale. The experimental data shown 
here were obtained also as a numerical derivative of the 
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Mea;/kJmol ^ A/^ei/kJ mol 



Model 


T/K 


Nc 


Ar 


Kr 


Xc 


Xct 


Xc* 


CH4 


Xc 


SPC 


275 


10.80 


7.64 


6.88 


6.18 




3.39 


8.14 


-2.59 




300 


11.41 


8.68 


8.12 


7.65 




4.68 


9.39 


-3.03 




325 


11.76 


9.32 


8.86 


8.44 




5.56 


10.06 


-3.61 




350 


11.96 


9.77 


9.38 


9.04 




6.22 


10.54 


-3.95 




375 


12.01 


10.07 


9.75 


9.45 




6.71 


10.84 


-4.35 


SPCE 


275 


10.95 


7.76 


7.06 


6.46 




3.59 


8.38 


-2.21 




300 


11.73 


8.97 


8.41 


7.99 




5.11 


9.71 


-2.82 




325 


12.31 


9.94 


9.54 


9.27 




6.08 


10.82 


-3.24 




350 


12.70 


10.61 


10.37 


10.24 




7.32 


11.58 


-3.71 






iz.yu 


ii.U ( 


iU.oO 


iU. 1 1 




( . ( 4 


iz.Uo 




SPCE (500 Mol.) 


275 


10.89 


7.66 


6.92 


6.15 


6.16 


3.36 


8.18 






300 


11.69 


8.87 


8.29 


7.71 


7.72 


5.00 


9.56 






325 


12.27 


9.84 


9.41 


9.05 


8.96 


6.22 


10.67 






350 


12.64 


10.52 


10.17 


9.87 


9.88 


7.03 


11.39 






375 


12.84 


10.94 


10.66 


10.41 


10.42 


7.56 


11.83 




SPCE (864 Mol.) 


275 


10.88 


7.62 


6.80 


6.01 




3.66 


8.09 






300 


11.65 


8.83 


8.20 


7.58 




4.87 


9.45 






325 


12.24 


9.75 


9.24 


8.74 




6.11 


10.54 






350 


12.61 


10.44 


10.05 


9.65 




6.91 


11.26 






o7o 


12.82 


10.91 


10.63 


10.34 




7.58 


11.78 




TIP3P 


275 


10.94 


7.85 


7.10 


6.40 




2.71 


8.34 


-2.81 




300 


11.39 


8.69 


8.14 


7.60 




4.09 


9.33 


-3.30 




325 


11.63 


9.25 


8.79 


8.35 




4.85 


9.97 


-3.74 




350 


11.72 


9.54 


9.11 


8.69 




5.37 


10.24 


-4.00 




375 


11.65 


9.66 


9.28 


8.88 




5.77 


10.34 


-4.32 


TIP4P 


275 


10.88 


7.77 


7.09 


6.40 




4.06 


8.38 


-2.52 




300 


11.64 


8.99 


8.49 


8.16 




5.82 


9.78 


-3.13 




325 


12.14 


9.81 


9.44 


9.28 




6.67 


10.71 


-3.54 




350 


12.39 


10.31 


9.97 


9.70 




7.17 


11.14 


-4.16 




375 


12.44 


10.59 


10.31 


10.06 




7.53 


11.44 


-4.43 


TIP5P 


275 


9.76 


6.36 


5.53 


4.79 




2.00 


6.81 


-1.36 




300 


10.79 


7.84 


7.19 


6.47 




3.34 


8.42 


-2.46 




325 


11.32 


8.76 


8.23 


7.76 




4.47 


9.44 


-3.25 




350 


11.45 


9.11 


8.62 


8.16 




5.15 


9.78 


-3.83 




375 


11.30 


9.16 


8.70 


8.25 




5.33 


9.78 


-4.02 



TABLE IV: Calculated excess chemical potentials fiex for the different solutes along the 0.1 MPa isobar obtained by the particle insertion 
technique. The Xe* column is obtained by taking the polarisability of o = 4.11 into account. The Xe^ column contains data accordig 
to the overlapping distribution method. The data in the most right column represent the change of the chemical potential when bringing 
a particle from infinity to the distance of the maximum of the Xe-Xe radial distribution function: A/iex = /iei: (0.42 nm) — fj,ex{oo). 
The accuracy of the data was estimated to: Ne : ±0.05 kj mol"^ ; Ar : ±0.08 kj mol~^ ; Kr : ±0.1kjmol~^; Xe : ±0.15 kj mol"'^; 
Xe* : ±0.25kJmori; Xet : ±0.05 kj mol"!; CH4 : ±0.1kJmol-\ A^lea:{Xe) : ±0.15kjmol-^ 



experimental excess chemical potentials with respect to 
the temperature. 

A detailed look at the experimental data for Methane 
according to Fernandez-Prini and Crovetto ■81] reveals 
decreased slope of s^^ in the region below 300 K. The 
reason for this is not quite clear and not present when cal- 



culating entropies calculated from the solubility data of 
Rettich et al. Since this observation fo^l is also not com- 
patible with the calorimetric measurements of Naghibi et 
al. it may probably reflect a deficiency of the fit used 
in Ref. jsj- Above 300 K both experimental datasets, 
however, do agree well. 
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The solvation entropy obtained for Methane in TIP4P 
water at 300K is found to be — 47±5 JK^^moP^, which 
is shghtly larger than the value of — 40.7 J K^^moP^ re- 
ported by Shimizu and Chan. The entropies are through- 
out negative, in accordance with the classical interpreta- 
tion of hydrophobic hydration, suggesting an enhanced 
ordering of the molecules in the solvation shell. In the 
region between 275 K and 300 K the experimental data 
reveal about 15— 10 JK^^moP^ larger absolute solvation 
entropies for Xenon compared to Methane. This trend is 
also observed in the simulations, where we find 7, 8, 5, 7 
and 4 JK~^moP^ larger absolute solvation entropies for 
Xenon for SPC, SPCE, TIP3P, TIP4P and TIP5P-water, 
respectively. However, the different models show a clear 
hierarchy with respect to their entropy when compar- 
ing with the experimental data for Methane and Xenon. 
For both cases. Methane and Xenon, TIP3P and SPC 
exhibit the smallest solvation entropies, whereas TIP4P 
and SPCE are lying closer to the experimental data. The 
TIP5P model, however, shows the strongest temperature 
dependence, exhibiting the about smallest solvation en- 
tropy at 275 K of all models and the highest value at 
375 K. The temperature dependence of the solvation en- 
tropies can be quantified in terms of solvation heat capac- 
ities. The decrease of the absolute value of the solvation 
entropies leads to positive solvation heat capacities (cal- 
culated for 300 K) which are estimated to be (145 ± 20), 
(144±20), (145±20), (184±20), (264±20) JE'VoP^ 
for Methane and (156 ± 20), (164 ± 20), (160 ± 20), 
(223 ± 20), (310 ± 20) JK-VoP^ for Xenon in SPC-, 
SPCE-, TIP3P-, TIP4P- and TIP5P-water, respectively. 
The experimental data of cp^^x according to the solubility 
data of Methane of Rettich et al. |63| was obtained to be 
234JK"^moP^ at 300 K, when transforming their data 
on the number density scale. This value has also been re- 
ported by Widom et al. 47] . The value of cp^ex for Xenon 
according to Fernandez-Prinis data using the same pro- 
cedure was estimated to be 280JK^^moP^. The gen- 
eral trend, larger absolute solvation entropies and heat- 
capacities for Xenon compared to Methane, is accom- 
plished by all water models and is also found in exper- 
iment. Moreover, the experimental data indicate an in- 
crease of the heat capacities at lower temperatures. The 
fitted curves shown in Figure |3| tend to suggest this as 
well. However, the error in the entropies obtained from 
finite differences is perhaps too large to definitely confirm 
this. The values for cp^ex given above reflect an average 
considering the whole temperature range. 

We would like to point out that the apparent hierarchy 
of the models with respect to their ability to reproduce 
waters structure (00-correlation) and the temperature 
of maximum density is also observed when comparing 
the solvation entropies of simple solutes. TIP3P has by 
far the smallest solvation entropy and the SPC model is 
only slightly better. SPCE is apparently the best sol- 
vent model at higher temperatures (from the point of 
view of solvation entropies), whereas the TIP5P model 
is is closer to the experimental data below 300 K. Ap- 
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FIG. 4: Contribution of the polarisability to the excess chemical 
potential of Xenon in water. 



parently the water models that exhibit a more strongly 
pronounced structure are also subject to an enhanced or- 
dering of the molecules in the solvation shell, which is re- 
flected by more negative hydration entropies. However, 
this is perhaps not too surprising, since a water model 
that provides a better representation of waters structure 
might as well lead to a more realistic structural descrip- 
tion of the hydrophobic hydration shell. The reason for 
this might be related to the observation that the liquid 
water structure contains cavities suitable for hydrophobic 
particles js^. This is, of course, also the cause for the 
Widom particle insertion technique performing so well 
for hydrophobic particles in water. 

The density-curves and corresponding structural fea- 
tures of the individual water models provoked the sim- 
plistic interpretation that TIP3P and SPC corresponds 
structurally to states of water at increased temperature. 
We would like to point out that the simulated entropies 
support this point of view. To a flrst order approxima- 
tion the simulated entropy-data can be simply corrected 
by a temperature shift. 

Finally, we would like to discuss the effect of hav- 
ing an explicitly polarisable solute particle. As shown 
in Figure 01 introducing a polarisable Xenon particle 
leads to a lowering of the chemical potential of about 
2.5kJmoP^ to 3.5kJmoP^. Please note that the con- 
tribution to the chemical potential due to the polaris- 
ability is only weakly temperature dependent. The most 
extreme cases in this respect are the TIP3P and the 
TIP4P model. The TIP3P-data show a change of about 
0.5kJmoP^ over the entire temperature range, indicat- 
ing a lowering of the entropy for each temperature on av- 
erage of about 5 JK^^moP^ For the TIP4P model, how- 
ever, this contribution is about — 2JK~^moP^, leading 
to a slight increase of the solvation entropy. Both values 
are lying in the range of the expected error for the en- 
tropy values derived from finite difference values of about 
±8 JK~^moP\ 
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C. A Modified Information Theory Model 

In the previous section we discussed the properties 
of hydrophobic hydration by relating structural features 
and solvation entropies employing (to a first order ap- 
proximation) a temperature shift. In this section we will 
try to elaborate a more quantitative way to describe the 
differences between the models. 

The information theory approach to the hydrophobic 
hydration according to Hummer et al. suggests that 
the excess chemical potential of a hard sphere particle 
(i.e. a spherical cavity) can be expressed as 



2 2 



(13) 



where v is the volume of the spherical cavity, p is the 
water number density and ct„ = ("'^) ~ (")^ is the vari- 
ance of the occupancy number of water molecules in the 
sphere. For the case of noticeable attractive interactions 
Hummer et al. [s^ l propose the presence of an additional 
attractive term plus an offset according to 



^ « Ap'-Cp/T + B 



(14) 



where A = v p''/{2a^) and the parameter C is found to 
qualitatively account for the effects of attractive solute- 
solvent interactions 83]. However, in order to arrive at 
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FIG. 5: Scaling plots of the excess chemical potential of Methane 
and Xenon in water according to the MIT-model. The MIT- model 
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data shown here, assuming linearity between fiex and Vm/T. 
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a quantitative description of the experimental and sim- 
ulated data, we find it advantageous to express the pa- 
rameters A and C, both inverse proportional to the tem- 
perature with 



A = a'/T and C c'/T . 



(15) 



The justification for this approach is purely empirical but 
can be perhaps rationalised as: 1.) An effect of a temper- 
ature dependent effective particle diameter, as well as a 
slight increase in the fluctuation cr^ 's^ . At higher tem- 
perature, the effective diameter of the solute is likely to 
decrease due to the form of the repulsive potential. 2.) 
An effect of a more strongly weakened attractive interac- 
tion with increasing temperature. Moreover, when doing 
this, the offset B can be dropped {B = 0), although we 
must admit that the original model in Eq. 1 141 using three 
parameters represents the data slightly better. 

Since this fitting procedure has been basically inspired 
by the information theory (IT) approach for hydrophobic 
hydration and interaction we will refer to it as modified 
information theory model (MIT) in the course of this pa- 
per. Moreover, the parameters a and c are expressed in 
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Model 


MIT-parameters 


Ne 


Ar 


Kr 


Xe 


Xe* 


CH4 


Xe" 


SPC 


a/lQ-^KinVor^ 


1.130 


1.169 


1.213 


1.254 


1.050 


1.272 


0.640 




c/iV m moi 


-lU. / 




1 /I Q 
-14. o 


-10. o 


-lo.y 


1/1/1 

-14.4 


1 A1 

- i Ai 


SPCE 


a/lO^'^KmVor^ 


1.085 


1.159 


1.218 


1.280 


1.076 


1.271 


0.749 




c/xv m mol 


-iU.o 


-io.4 


-14. / 


1 f; A 
-ib.U 


-14.0 


-14. ( 


AQ 

-y.uo 


TIP3P 


a/10"^KmVol-2 


1.121 


1.138 


1.165 


1.184 


1.002 


1.230 


0.569 




c/iv m moi 


-iU.o 


1 o o 


-io. / 


1/1/1 

-14.4 


-lo.o 


1 /I n 
-14. u 


-0.04 


TIP4P 


a/10"''Km^mor^ 


1.096 


1.151 


1.190 


1.240 


1.029 


1.245 


0.654 




c/K m mol 


1 A C 

-iU.o 


1 O 1 

-13. i 


1 yl 1 

-14.1 


ICO 

-15.2 


-13.3 


1 /I 
-14.2 


-7.61 


TIP5P 


a/10-'^KmVol-2 


1.133 


1.173 


1.200 


1.223 


0.972 


1.256 


0.475 




c/K^m^mol"^ 


-11.3 


-13.9 


-14.8 


-15.6 


-13.5 


-14.9 


-5.08 


Expt. 


a/lO^^KmVor^ 


1.104 


1.191 


1.212 


1.174 




1.243 






c/K^m^mol-i 


-10.9 


-14.2 


-15.5 


-15.7 




-15.1 





TABLE V: Parameters describing the excess chemical potential of the noble gases and Methane applying the modified information theory 
(MIT) model. The model parameters were obtained by fitting the data of Table ITVl The Xe*-data corresponds to the polarizable Xe- 
particle. The Xe'^-data corresponds a Xenon par ticl e in c onta ct with another Xenon particle. The parameters representing the experimental 
data were fitted to the data published in Refs. l8ll and l7ll . 



terms of the molar volume, such that 



with R being the ideal gas constant. 

In Figure |5l we show scaled plots of the chemical po- 
tential of Xenon and Methane according to eauationll6l 
Taking the experimental data, the two parameter MIT 
model accurately represents the data of the noble gases 
up to temperatures of 420K — 450K. At higher temper- 
atures the linearity between /iea; and Vm/T breaks 
down, which is not unexpected due to the enhanced in- 
crease of the isothermal compressibility kt ^3- Com- 
paring the obtained MIT-parameters for the noble gases 
(see Table W\ . it is evident that the parameters also be- 
have meaningfully when going from Ne to Xe, in the sense 
that the c-parameters becomes more negative, hence the 
attractive interaction becomes stronger, whereas the in- 
creasing a-parameter accounts for the increasing size of 
the hydrophobic particle. 

Nevertheless, the limits of the model become evident 
also: The slight deviation from linearity for Methane 
(which is fully absent in the case of Xenon) in Figure 
El becomes significantly stronger in the case of Neon (not 
shown), indicating that both parameters are subject to 
a enthalpy-entropy compensation effect and a sufficiently 
strong attractive (e.g. a large Lennard- Jones e) interac- 
tion is required for a good performance of the model. Or 
in other words: A counterbalancing of the two terms is 
necessary. 

Moreover, it is observed that the a parameter obtained 



for Xenon is even smaller (or at least very close) to the 
value obtained for Methane. A stronger attractive inter- 
action apparently leads to a shrinking apparent cavity 
size which is maybe related to the presence of a more 
tightly bound hydration shell. This feature is strongly 
pronounced in the case of the experimental values, here 
the a-parameter for Xenon is even smaller than the value 
for Krypton. This might be well attributed to the larger 
polarisability of the Xenon atom compared with Krypton 
and Methane. As shown in Table [V] taking the polaris- 
ability into account significantly reduces the a and c. 

Figure|Slreveals that the scaled chemical potential data 
for all models (except for TIP5P) arc lying quite close to 
each other. In addition, the experimental line has almost 
the same slope, so that the difference between the dif- 
ferent models (and experiment) can be described almost 
completely by just shifting the a-parameter. Moreover, 
the data according to the different models (except for 
TIP5P) fall almost onto one line. As a consequence it is 
worthwhile trying to describe the chemical potentials by 
using the same set of parameters (here we take the data 
corresponding to the SPCE-model), while just taking the 
density-curves of the different models into account. 

In Figure Et- the Xenon excess chemical potentials for 
the different water-models, as well as the MIT-model 
predictions based on the SPCE-parameters are shown. 
An almost quantitative prediction is achieved, except 
for the TIP5P model. A possible explanation for the 
large deviation of TIP5P-data might be the noticeably 
smaller Lennard- Jones a of the TIP5P-model (see table 
P) . The information theory suggests that the a-parameter 
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should scale with the square of the the particle volume. 
Hence we try to improve the model prediction by scaling 
the a-parameter with the factor (crxe-modoi/cxc-SPCE)^, 
shown in Figure]^, which, indeed, leads to a substantial 
improvement for the case of TIP5P. In order reproduce 
the experimental data, a scaling procedure taking the 
experimental density-curve and and Lennard-Jones pa- 
rameter of (Txc-Watcr = 3.5275 A has to be employed. 

Apparently the suggested MIT model is able to de- 
scribe the excess chemical potential of Xennon for the 
different water models and experiment by just taking 
the different density isobars into account. In the previous 
section, however, it was argued that both water structure 
and the solvation entropies are well related and hence 
responsible for the performance of the different models. 
However, this is perhaps not a contradiction when keep- 
ing in mind that waters structural and density changes 
are tightly related and that the transformation towards 
a more tetrahedrally ordered structure at lower tempera- 
tures is the basis for the presence of a density maximum 
0- 

We are quite confident that similarities between ex- 
perimental and simulated data applying the suggested 
rescaling procedure reveal a signature for hydrophobic 
hydration in the lower temperature regime and might 
also be helpful when trying to quantify the interaction 
between hydrophobic of particles. 
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D. Hydrophobic Interaction 

In order to quantify the hydrophobic interaction we 
calculate the Xenon- Xenon pair distribution functions 
for the different models and temperatures. The data 
are shown Figure |H1 AH models show an increase of 
the first peak with rising temperature. This observa- 
tion is well in accordance with the interpretation that 
the association of two hydrophobic particles is stabilised 
by minimising the solvation entropy penalty and has al- 
ready been reported by a large n umber of publications 
Hllllllllliailllllll I2II23 ■ The idea is schemat- 
ically depicted in Figure [7| The enhanced ordering of 
the solvent molecules in the hydration shell results in 
a negative solvation entropy. If two particles associate, 
the corresponding hydration shells overlap, hence lead- 
ing to a positive net entropy. Consequently contact- 
configurations should become increasingly stabilised at 
higher temperatures. In parallel, the increased heat ca- 

( ; o) ^ (00) 

^ 

FIG. 7: Schematic diagram of the hydrophobic association pro- 
cess according to the classic picture: The contact configuration is 
stabihsed with increasing temperatures by minimising the entropy 
penalty. 



FIG. 8: Xe-Xe radial pair distribution functions g{r) for the dif- 
ferent water models and temperatures. The arrow indicates the 
sequence of g(r)-curves pointing from low to high temperatures. 
The bar indicates the change of the height of the first maximum 
over the whole temperature range. 



pacity of the water molecules in the hydrophobic hydra- 
tion shell should lead to a negative heat capacity con- 
tribution for the association of two particles, thus weak- 
ening the the entropy contribution at elevated temper- 
atures. We would like to emphasise that this model 
of course lacks completely detailed structural consider- 
ations. Our study reveals significant differences for the 
different water models. The overall maximum variation 
of the height of the first peak (which is also indicated 
in Figure IHl) is smallest for the TIP3P and SPC models, 
TIP4P and SPCE show a stronger variation, whereas the 
TIP5P model reveals an extremely pronounced dropping 
of the first peak when going to lower temperatures. In 
all cases the maximum of the first peak is located at a 
distance of about 4.2 A, which does not change signifi- 
cantly with varying temperature. Moreover, the pair dis- 
tributions functions reveal as well the presence of a pro- 
nounced second peak corresponding to the solvent sepa- 
rated Xenon-Xenon pair configuration, which is located 
at 7.2 — 7.8 A and is shifting to smaller distances with 
decreasing temperatures. In addition, a third maximum, 
being located at about 10.5 A, seems to exist at least at 
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FIG. 9: Cavity part of the profile of free energy w{r) obtained 
for the association of two Xenon particles for three selected models 
and all temperatures. The heavy dashed lines denote the Xe-Xe 
Lennard-Jones pair potential. 




FIG. 10: Entropy profile s{r) for the association of two Xenon 
particles. 



the lowest temperatures. 

In Figureinithe corresponding profiles of free energy for 
the association of two Xenon particles are shown. The 
TIP3P-, SPCE- and TlP5P-water models might be taken 
here as the representative cases for small, medium and 
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FIG. 11: Enthalpic and entropic contributions to the cavity part 
of the profile of free energy w{r) = h{r) — Ts{r). Solid lines: h{r). 
Dotted hues: -Ts{r). 



strong temperature variation, respectively. In addition, 
the relative change of the excess chemical potentials when 
bringing a Xenon particle from the bulk to the Xenon- 
Xenon contact distance A/^tej; = ^,f,x{QA2mn) — ^f,x{oo) 
are given in Table IIVI for all models and temperatures. 
We would like to point out that quadratic fits with re- 
spect to the temperature describe the changes of data 
over the entire temperature range reasonably well. From 
the distribution of pair correlations functions obtained 
from different parts of the simulation runs, we conclude 
that the individual curves shown here are accurate within 
an interval of ±0.15 kJ moP^. For completeness Figure|51 
contains as well the temperature independent pair poten- 
tial between two Xenon particles. Figure |51 shows that 
with decreasing temperature the contact configuration is 
increasingly destabilised, whereas the solvent separated 
configuration becomes more and more stable. 

From the temperature derivative of the w(r, T) data 
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set we obtain the entropy-profiles for the association of 
two Xenon particles, shown in Figure IIUI The entropic 
—Ts{r) and enthalpic h{r) contributions to the profile 
of free energy are given in Figure 1111 The temperature 
variation of the entropy profiles is quantified by the corre- 
sponding heat capacity profiles cp(r) given in Figure [HI 
As already shown by Smith and Haymet 13, P| and others, 
the hydrophobic association process is found to be en- 
tropically favoured and enthalpically disfavoured in case 
of all models. The value of —5.3 kJ mol~^ for the entropic 
contribution to the profile of free energy at 300 K at the 
contact distance for TIP3P water is reasonably close to 
the « — 3.5kJmol~^ (when taking the data from Fi gure 
5 in Ref. [23|) reported for Methane by Ghosh et al. |27l| 
and the value of — G.SkJmoP^ is in the same sense con- 
sistent with the -4.14kJmor^ for Methane in TIP4P 
water obtained by Shimizu et al.|2^. The observation 
of slightly larger entropies for the association of Xenon 
instead of Methane is plausible since the hydration en- 
tropies for Xenon are also found to be more negative. We 
would also like to point out that the hierarchy of the data 
for the different models according the published data is 
apparently consistent with our observations. The data of 
Rick (24l | obtained from simulations employin g th e polar- 
isable TIP4P-FQ model have been criticised |23ll2l l27t 
because of the unreasonably high heat capacity change 
for the association of two Methane particles. However, 
the large entropy contribution to the profile of free energy 
of about — 12kJmoP^ at the contact distance is rather 
close to the value of — ll.SkJmoP^ that is observed here 
for the TIP5P model. 

Figure EH shows that in all cases there is a ten- 
dency for the entropy at very short distances to become 
smaller with increasing temperature, whereas in the re- 
gion around 6 A, at the so called desolvation harrier, the 
entropy increases with temperature. The general ten- 
dency to smaller contact-entropies is in accordance with 
the decrease of the absolute values of the solvation en- 
tropies. However, the distances where the entropy-curves 
are found to cross each other are quite diverse. Whereas 
for the TIP-models this region is located around 5.5 A it 
is found in case of the SPC and SPCE models at about 
4.5 A, much closer to the Xenon-Xenon pair distance of 
about 4.2 A. As a consequence the pair formation en- 
tropies vary much more strongly in case of the TIP mod- 
els. 

From Figure ^] it is evident that in almost all cases 
the pair configuration is entropically stabilised and en- 
thalpically destabilised. However, TIP5P shows a much 
stronger enthalpy-entropy compensation effect than all 
other models, but as well a much stronger variation with 
temperature. In addition, also the signature of the hy- 
drophobic association vanishes in case of the TIP5P and 
TIP3P models at 375 K where the entropic and enthalpic 
contribution to the cavity potential at the the contact 
distance have almost the same size. This is consistent 
with finding that in case of the TIP5P and TIP3P mod- 
els maximum of ^ex for Xenon is close to 375 K and with 



the contact entropy Sex being close to zero. 

The temperature dependence discussed here can be 
quantified by the association heat capacities cp(r) shown 
in FigurelT^ Shimizu and Chan 23, 25, 85] report for the 
association of Methane particles in TIP4P water a change 
of the heat capacity for the contact state close to zero in 
qualitative disagreement with the approximate net effect 
according to the overlapping hydration shells. In addi- 
tion, they observe a maximum of the heat capacity of 
about 120 J mol~^ at the location of the desolvation 
barrier at a distance of 5.5 A. In a previous study using 
a polarisable water model Rick |0| did not observe any 
peak at the desolvation barrier. Nevertheless, his value of 
about —2500 J mol~^ for the relative change of heat 
capacity for the contact state is unreasonably large, being 
about more than one order of magnitude larger than the 
solvation heat capacity of Methane found in experiment 
and the water models used in our study. However, in a 
more recent study Rick ^86j confirms the presence of a 
peak at the desolvation barrier of about 40 J K^^ mol~^ 
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FIG. 12: Relative change in lieat capcacity cp(r) for the hydropho- 
bic interaction between two dissolved Xenon particles. The data 
correspond to the temperature derivative of quadratic fits of u)(r, T) 
obtained for 300 K. For displaying purposes the data are shifted 
by an offset of 50 J mol~^. 
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a peak in the heat capacity at the desolvation barrier in 
the region around 6 A of about 20 JK~^niol~^ for the 
TIP3P and SPC model, about 70 JK^VioP^ for the 
TIP5P and SPCE model and about 50 JK^VoP^ for 
the TIP4P. 

However, for the change in heat capacity for the con- 
tact state the different models show a quite diverse be- 
haviour: All TIP models reveal a negative contribu- 
tion to the heat capacity for the contact state (TIP3P: 
-44JK-^mor^; TIP4P: -42JK-^mor^; TIP5P: 
-130 JK"^ mol"^) with the value for the TIP5P model 
being extremely large. The Berendsen models reveal sig- 
nificantly smaller values with —23 JK~^mol~^ for SPC 
and -1-5 J mol~^ for the SPCE model. Hence here we 
apparently observe for Xenon and the SPCE model what 
Shimizu and Chan observe for Methane in TIP4P water: 
A slightly positive heat capacity contribution. The ap- 
parent disagreement between the behaviour of Methane 
and Xenon in TIP4P water, however, might be related to 
the detailed hydrogen bonding situation (and their flucu- 
ations) around the differently sized particles. We would 
also like to point out that Rick 86] finds a negative heat 
capacity contribution for the Methane-Methane contact 
pair in TIP4P water of about -140 ± 120 JK"^ moP^ 

The apparently smaller heat capacity contributions 
of the Berendsen-models and the extremely large value 
found for the TIP5P model, however, strongly suggest 
that details in the hydrogen bonding situation (i.e. struc- 
ture of the hydration shell) of the Xenon-Xenon contact 
state are responsible for the differences. A further de- 
tailed comparative analysis concerning the energetics, hy- 
drogen bonding fluctuations and structure of the hydra- 
tion shell for the different models will help to reveal this 
differences more clearly and is the topic of a following 
publication. 



repulsive water cavity potential. In addition, the desta- 
bilisation of the contact-state with decreasing tempera- 
ture is even enhanced by the increased lowering of the 
solvent separated minimum of the profile of free energy 
as shown in Figure |51 The calculated intersection tem- 
perature for TIP5P water is found to lie quite high with 
257 K. For the other water models, however, the predic- 
tion of the intersection temperature is much more prob- 
lematic since the temperatures are apparently shifted to 
even lower values. Hence the density fits, necessary to 
determine the intersection temperatures, are much less 
precise. Therefore the values of about 220 K for SPC 
and TIP3P and 230 K for SPCE and TIP4P arc subject 
to large uncertainty. 

In section UlI CI we have shown that the experimental 
data as well as the TIP5P and SPCE data can be almost 
quantitatively interrelated by the MIT model assuming 
an effective Lennard Jones sigma for the experimental 
water of crxe-Water = 3.5375 A. When applying the same 
procedure here, however, we should get an impression of 
how the strength of the hydrophobic interaction of Xenon 
in real water might behave. We would like to point out 
that the SPCE and the TIP5P model represent the most 
extreme cases with respect to their slopes shown in Fig- 
ure El The so calculated strengths of the hydrophobic 
interaction as a function of temperature are shown in Fig- 
ure ll4l and compared with the corresponding data for the 
water models. Besides the large uncertainties of this ap- 
proach, the most striking feature, the strongly changing 
slope of the A/Xe^, curve, is mostly due to the tempera- 
ture dependence of waters expansivity. Hence this obser- 
vation does not depend on the exact values of the a and c 
parameters, but requires just an approximately linear de- 
pendence between fi^x Kn ^'^'^ Vm/T, as it is suggested by 
all five water models. The corresponding temperatures 



E. Concerning the Hydrophobic Interaction 
Between Xenon Particles in Real Water 



In Figure [T^ we have applied the rescaling procedure 
proposed in section IIII CI to the chemical potentials of 
the Xenon particles in contact with another Xenon par- 
ticle. Again we obtain an apparently linear behaviour 
for each of the water models, although we must admit 
that the corresponding slopes are found to vary more 
strongly than for the case of the bulk liquid. However, it 
is quite evident that all slopes are considerably smaller 
than the ones obtained for the bulk. The fit parameters 
are given in Table IVl An obvious conclusion is of course 
that the lines corresponding to bulk and shell might have 
an intersection, defining the temperature where the in- 
teraction between the hydrophobic particles turns from 
attractive into repulsive. The strongest evidence for this 
scenario comes from the TIP5P data with a value of 
A/iex(0.42 nm) = — 1.36kJmol~^ at 275 K, which is ac- 
tually weaker than the pure Lennard Jones attraction 
between the two Xenon particles, indicating an already 
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FIG. 13: Scaled plots of the excess chemical potential fiex of 
Xenon in water for the different water models according to the 
MIT-model. Open symbols: chemical potential obtained for the 
water bulk. Closed symbols: chemical potential in the vicinity of 
a Xenon particle (obtained at the distance of the first peak of the 
Xe-Xe g(r)-function with r= 0.42 nm). The lines represent MIT 
fits taking all water models into account. The data corresponding 
to the lowest temperatures are found on right side of this diagram. 
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FIG. 14: Excess chemical potential of Xenon in the vicinity of an- 
other Xenon particle relative to the bulk. The distance corresponds 
to the maximum of the Xe-Xe 3{r)-function (r = 0.42 nm). The 
heavy dashed and dot-dashed lines represents the model prediction 
for the association of two Xenon particles in real water based on 
the MIT-parameters for SPCE and TIP5P (see text for details). 
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FIG. 15: Variation of the strength of the hydrophobic interaction 
A/i(0.42nm) = fj,ex{0A2nm) — fiexioo) with the hydration free en- 
ergy flex- The heavy dashed and dot-dashed lines represents the 
model prediction for the association of two Xenon particles in real 
water (shown in Figure lMl based on the MIT-parameters for SPCE 
and TIP5P. 



for the attractive/repulsive conversion are hence found 
to He in the interval between 253 K and 267 K. 

Recently Widom et al. [i^l have deduced from lat- 
tice model calculations the presence of an almost linear 
relation between the hydrophobic interaction Aii^^/RT 
and the free energy of hydrophobic hydration ^^x/RT 
for the limited temperature interval between 273 K and 
333 K. Figure [T51 shows a such a plot containing the data 
obtained here for the different water models. Figure ITKl 
indicates that at least for temperatures sufficiently below 
the maximum of (T) / RT their observation is consis- 
tent with our data. For the case of the SPCE, TIP4P and 
TIP5P models an almost linear behaviour in the interval 
275K and 325K is observed. Please note that the model 
predictions for "real Xenon in water" based MIT model 
as outlined above behave as well approximately linear for 
the lower temperatures, suggesting that the MIT predic- 
tions for the lower temperatures are as well conceptually 



consistent with the theory of Widom et al. [Tg]. We 
would like to point out that the slope of 0.7 observed by 
Widom et al. for Methane is closer to the value of 0.8 
obtained for the TIP5P model than to the « 0.4 found 
for SPCE and the other water models. 



IV. CONCLUSIONS 

From a series of Molecular Dynamics simulations on 
different water models we conclude that the differences 
between the model waters structure, expansivity behav- 
ior and the temperature dependence of the solubility of 
simple solutes are tightly related. The calculated solva- 
tion entropies for all models are found to be systemati- 
cally smaller than the corresponding experimental values. 
However, the water models (SPCE, TIP4P, TIP5P) that 
reproduce waters structure more accurately, provide sol- 
vation entropies that are closer to the experimental data. 

According to a modification of the Information theory 
model of Hummer et al. we observe an almost linear de- 
pendence between n^x ^'^'^ Vm/T for both experimen- 
tal and calculated excess chemical potentials. A corre- 
sponding rescaling procedure is able to almost eliminate 
the density dependencies of the different data. Moreover, 
when taking the size of Xenon particle according to the 
different Lennard-Jones sigmas for the different models 
into account, an almost quantitative prediction of the ex- 
cess chemical potential of Xenon for the different models 
and the experimental data is feasible. 

Concerning the hydrophobic interaction we would like 
to emphasise the simple fact that all models show a qual- 
itatively similar behaviour: Enhanced aggregation at ele- 
vated temperatures. This is in general consistent with the 
simple picture of the association being driven by entropy 
effects determined as at net result of overlapping hydra- 
tions shells. Nevertheless, from a quantitative point of 
view noticeable differences between the different models 
exist, which can only be rationalised considering molec- 
ular details of the hydration shell in the Xenon-Xenon 
contact state. However, in general we can state that the 
water models that reveal a more realistic hydration en- 
tropy and water structure, namely SPCE, TIP4P and 
TIP5P, reveal as well a more strongly pronounced associ- 
ation behaviour. The TIP5P model reveals an extremely 
pronounced temperature dependence. 

The apparent linear behaviour between fi^x and 
Vm/T is also found for Xenon in the Xenon-Xenon con- 
tact state, strongly suggesting the existence of a tem- 
perature where the hydrophobic interaction turns from 
attractive into purely repulsive. For the case of TIP5P 
water this temperature is found to lie quite high at about 
257 K. Moreover, the strong temperature dependence of 
(experimental) waters expansivity strongly indicates that 
the weakening of the hydrophobic interactions in real 
water is closer to that obtained for the TIP5P model. 
Consequently, a model that accounts for waters density 
effects as closest as possible is as well desirable for a cor- 
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rect description of hydrophobic interactions. 

The almost linear relationship between hydrophobic 
interaction and the stren gth of hydrophobic hydration 
proposed by Widom et al. M7| is at least for temperatures 
sufhciently below the maximum of (T) / RT consistent 
with our data. 
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APPENDIX A: APPENDIX: CHECK ON THE 
APPLICABILITY OF THE PARTICLE 
INSERTION TECHNIQUE BY THE 
OVERLAPPING DISTRIBUTION METHOD IN 
THE ISOBARIC-ISOTHERMAL ENSEMBLE 

The Widom particle insertion method is known to fail 
in some cases, such as e.g. the high density Lennard- 
Jones liquid |6j]. In order to confirm that the Widom 
method can be applied under the conditions of our study 
we have used additional simulations of one Xenon particle 
in 500 water molecules to calculate the chemical poten- 
tial for Xenon using the overlapping distribution method 
as outlined in |6lL l88j . The additional simulations were 
conducted under the same temperatures/pressures and 
simulation parameters as outlined in section lll Al and ex- 
tended over 10 ns. Since Xenon is the largest particle in 
our study, it has the smallest free volume and therefore 
represents the worst case scenario comparing with other 
(smaller) noble gases and Methane. 

Since we are dealing with simulations in the isobaric- 
isothermal ensemble, the definition of the distribution 
functions has to be slightly different than that for the 
canonical ensemble, although this difference is usually 
ignored [s^ . In analogy to the formulation of the over- 
lapping distribution functions in the canonical ensem- 
ble |6l| we define two distribution functions (or his- 
tograms) pi and po of the energy of the -I- I'th particle 
AJ7 = C/(s^+i;L) - ;7(s^;L). Herepi(AJ7) denotes the 
distribution of energies of the A'' -I- I'th particle in the 
N + 1-particle system in the isobar isothermal ensemble 
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N+l j ua V X 
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exp(-/3Py) exp [-/3[/(s"+i;L)] x 
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scaled coordinates of particle N . and s^+^ repre- 
sent the set of coordinates of the entire iV-and -I- 1- 
particle systems, P is the pressure and V the volume. 
We have to denote that the partition function contains 
the factor A^! instead of (A^-|- 1)! here, since the A^+ I'th 
particle, the solute particle, is distinguishable from the 
N solvent particles in case it corresponds to a different 
particle type. U{s^;L) and U{s^^^; L) denote the po- 
tential energies of the A'' and A^ -|- 1 particle systems, 
respectively. Q{. . .) is the isobaric isothermal partition 
function. Eq. lAll is a straightforward extension of the 
expression for the canonical ensemble |6lj| . The func- 
tion po{AU) describes the energy-distribution of an ad- 
ditional A'^ + 1-particle randomly inserted into configu- 
rations representing an isobaric-isothermal ensemble of 
the A^-particle system. In contrast to the pi(AC/) dis- 
tribution function, the pq{AU) distribution, however, is 
defined with the instantaneous volume being used as a 
weighting factor 
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which has consequently to be normalised by the average 
volume {V). Starting with the definition of the pi(At/)- 
distribution function and inserting the definition for AC/ 
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where /3 = 1/kT and sn = L'^ tn (with L = V'^/^ be- 
ing the length of a hypothetical cubic box) represent the 
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FIG. 16: Excess chemical potential /lei according to the Widom 
particle insertion method and the method of overlapping distri- 
bution functions. All data were obtained at p = 0.1 MPa under 
NPT-conditions. The figure body shows /o, /ii /l — /o functions 
and the Widom estimate at T = 300 K. Insert: Well overlapping 
parts of /i — /o and Widom estimates for all temperatures. All data 
were obtained from simulations containing 500 SPCE molecules. 
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we come up with the following relation between the two 
distribution functions 



Pi(AC/) 



Q{N,P,T) (y) 
Q(N + 1,P,T) A3 
cxp(-/3AC/) po(AC/) 



(A3) 



Using the defition of the ideal and excess part of the 
chemical potential ^ referring to the ideal gas state with 
the same average number density in Eq. |21 we obtain a 
relation between the two distribution functions and the 
excess chemical potential which is analogous to the ex- 
pression for the canonical ensemble 

lnpi(A[/) - lnpo(AC/) = |3^ie^ - (3AU . (A4) 

We have to emphasise that the difference, however, is 
the necessity of volume-weighting in the calculation of 
the po(AJ7)-distribution function. As it is usually done 
|6l|. we define functions /o and /i using 



/o(a;7) = 
/i(a;7) = 

such that 



lnpo(AC/) - 
lnpi(AJ7) + 



f3AU 
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(3AU 



and 



h{AU)-fo{AU) 



(A5) 



In most cases the relative volume-fluctuation will be 
small, and therefore the volume- weighting will cause only 



a minor modification of to the po{AU) distribution. How- 
ever, when considering states close to the critical point, 
where volume fluctuations might become large, a signif- 
icant influence cannot be ruled out. Nevertheless, the 
above formulation should be used in any case, since the 
volume weighting can be done with practically no addi- 
tional computational effort. 

In Figure El the /i and /o distribution functions ob- 
tained from the 300 K simulation are shown, as well as 
difference between both functions and the excess chemi- 
cal potential obtained from the particle insertion method. 
The /i and fo distribution functions overlap largely and 
therefore a rather precise estimation of the excess chemi- 
cal potential of about ±0.05 kJmol"^ is feasible. Figure 
[TBI (see insert) shows a detailed comparison of /i — /q- 
data values according to the Widom method for all tem- 
peratures. The values obtained from the Widom particle 
insertion method agrees quantitatively with the data ob- 
tained from to the overlapping distribution method (see 
also Table Hvl for the data) The nice agreement between 
both methods, can be explained by the form of the /q- 
function. The /o exhibits a maximum, which indicates 
that the insertion procedure also explores the low-energy 
edge of the energy distribution with appropriate statis- 
tics. As a consequence there is no bias in sampling the 
energy space and the Widom method provides reliable re- 
sults. Employing the method of Flyvbjerg and Petersen 
j57j | we estimate an error of ±0.1 kJmol"^ as an upper 
bound for the accuracy of the Widom data for Xenon 
obtained from the 500 molecule system. 
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